/*--

 Based on an implementation of the algorithm MB-ID described in
 "Finding Exact Solutions to the Bandwidth Minimization Problem"
 by G. Del Corso  and G. Manzini
 modified by Jinwoong Hwang.

 I referred the source from the Corso and Manzini's implementation.
 However, I did not copy the whole code simply, but I reimplemented the code
 reading the paper.
 I assure that I understood the algorithm of the Iterative Deepening
 using expiration time concept.

 I adjusted some implementations and removed trivial codes
 and removed the pruning steps that I do not understand.


 This algorithm uses 'expiration' concept. 
 
 The idea starts from observation to the relationship between vertex degree
 and minimum bandwidth. The bandwidth cannot be less than 
 upper(max vertex degree /2) because it is linear - we can put adj nodes 
 on the left and the right.
 
 Iteration 1. 
 This iteration starts from b = upper(k/2) where k is max degree of the graph.
 and if b fails, it will try b+1, b+2 and so on. 
 
 Iteration 2. 
 Inside of the iteration 1, the algorithm tries adding nodes one by one 
 to the permutation and updates expiration time of nodes using DFS.
 In the procedure, there are two types of nodes. 
 One is free(busy) and the other is nonfree nodes. 
 When a node 'n' is selected, n's adjacent nodes will be 'hold' 
 as busy nodes. Also, every adj nodes will get decrement of their exp value.
 
 The queue of nodes keep being ordering by those exp value. So, adjacent nodes
 to the already selected nodes will be first in the queue. 
 By adding adjacent nodes one by one, a node close to the selected nodes will 
 be considered first, and have higher probability to get into the permutation.
 
 At each time after node choice is done, there should be checking compatibility
 process. If the choice is figured out that is not going to make it, the 
 procedure will backtrack to previous step. 
 
 If a candidate has expiration value 0, that means the candidate should be
 chosen in the step. If the thing that exp value falls below 0 happened, 
 it is obvious for the permutation to be failed.
 


SAMPLE RESULTS

//////// g-gg 
 
jin@ubuntu:~/Documents/bandwidth-minimizer$ ./band g-gg-36-60 
nodes: 36
edges: 60

Optimal ordering:
{ 1 2 27 36 15 14 25 33 22 28 10 29 11 13 19 16 35 34 18 5 24 32 9 30 21 26 31 8 6 12 17 20 7 23 4 3}
Minimum band: 6
Time:    33.7610(ms)

jin@ubuntu:~/Documents/bandwidth-minimizer$ ./band g-gg-91-162 
nodes: 91
edges: 162

Optimal ordering:
{ 5 23 26 29 8 48 6 31 39 18 68 34 1 84 58 81 66 62 67 79 53 86 75 37 69 35 56 51 49 3 88 59 12 71 21 38 90 45 73 30 10 70 15 47 13 19 41 50 32 46 36 2 55 54 60 63 85 44 16 61 72 78 9 14 82 4 77 7 57 43 52 42 27 28 24 20 33 11 65 87 76 64 91 25 83 40 89 17 74 80 22}
Minimum band: 7
Time:  9339.9560(ms)


//////// g-r 

jin@ubuntu:~/Documents/bandwidth-minimizer$ ./band g-r-30-109 
nodes: 30
edges: 109

Optimal ordering:
{ 3 9 21 4 19 30 29 8 17 11 2 25 26 12 24 20 7 16 10 22 1 18 6 14 28 15 5 27 23 13}
Minimum band: 13
Time:   718.1040(ms)

jin@ubuntu:~/Documents/bandwidth-minimizer$ ./band g-r-50-301 
nodes: 50
edges: 301

Optimal ordering:
{ 1 5 32 40 33 49 11 28 39 12 25 35 15 24 38 31 9 14 16 13 41 44 46 36 34 18 4 48 42 29 7 22 8 10 6 26 23 30 3 43 17 21 45 20 50 2 27 37 47 19}
Minimum band: 27
Time: 217916.9310(ms)


//////// g-bt

jin@ubuntu:~/Documents/bandwidth-minimizer$ ./band g-bt-20-19 
nodes: 20
edges: 19

Optimal ordering:
{ 12 6 13 15 3 7 14 1 5 10 2 11 20 4 9 18 8 19 16 17}
Minimum band: 3
Time:     7.1300(ms)


//////// g-p
jin@ubuntu:~/Documents/bandwidth-minimizer$ ./band g-p-100-99
nodes: 100
edges: 99

Optimal ordering:
{ 40 22 71 6 92 35 29 28 9 27 10 25 19 59 94 4 14 13 18 53 85 66 99 78 51 17 45 21 16 86 84 98 5 39 82 75 96 49 23 50 65 54 76 26 60 42 47 36 57 62 34 74 88 46 97 79 90 12 37 64 69 7 55 63 91 80 61 11 95 31 83 43 30 73 58 1 2 56 38 89 81 93 24 100 20 67 8 41 3 44 15 33 32 48 72 52 70 77 68 87}
Minimum band: 1
Time:     0.3250(ms)
jin@ubuntu:~/Documents/bandwidth-minimizer$ ./band g-t-40-39 
nodes: 40
edges: 39

Optimal ordering:
{ 2 10 38 24 15 32 33 37 3 18 13 28 17 23 40 4 29 5 12 21 16 1 39 11 8 35 27 36 30 34 9 20 7 6 19 31 26 14 22 25}
Minimum band: 3
Time:     0.1970(ms)


//////// g-t
jin@ubuntu:~/Documents/bandwidth-minimizer$ ./band g-t-40-39 
nodes: 40
edges: 39

Optimal ordering:
{ 2 10 38 24 15 32 33 37 3 18 13 28 17 23 40 4 29 5 12 21 16 1 39 11 8 35 27 36 30 34 9 20 7 6 19 31 26 14 22 25}
Minimum band: 3
Time:     0.1970(ms)


 
 --*/

#include <stdio.h>
#include <stdlib.h>
#include <time.h>  // clock_t, clock, CLOCKS_PER_SEC
#include <math.h>
//#include <sys/times.h> // POSIX ONLY
#include <sys/time.h>
#include <string.h>
#include <limits.h>
#include "bit.h"

#define PRINT_MATRIX 0
#define USE_TIME_LIMIT 0           // 1 to put a time limit
#define TIME_LIMIT_AMOUNT 120000   // time limit in milliseconds

#define MAX_SIZE 256    // Maximum size of martrix.
//#ifndef CLK_TCK 	// time.h code - something seems not to work in there.
#define CLK_TCK      CLOCKS_PER_SEC
//#endif

typedef short int Node;

typedef short int Nstep; 

/* ===================================================================
 For bit manipulation - it helps for performance
 =================================================================== */
#define lsize (31+MAX_SIZE)/32 // long to MAX_SIZE elements
typedef struct {
	unsigned long bits[lsize];
} list; // type used to represent a list. a row of adj list.

/* ===================================================================
 structure representing the current status
 =================================================================== */
typedef struct {
	Node queue[MAX_SIZE];	// list of nodes to be placed
	int qsize; 		// number of the nodes yet to be square
	int busy;		// number of nodes busy
	list notfree;		// compact list of nodes which are notfree
	Nstep expire[MAX_SIZE]; // list of expiration dates
} status;

/* ===================================================================
 structure containing the result of a test
 =================================================================== */
typedef struct result {
	double node_gen; // number of nodes generated
	double elapsed_time;
	int band; // bandwidth
} result;

/* =====================================================================
 matco - matrix compact - array of struct 'list' - fast bit manipulation
 matex - matrix expanded - char 2 dimension matrix
 ===================================================================== */
list matco[MAX_SIZE]; // representation of the matrix in compact form 
char matex[MAX_SIZE][MAX_SIZE]; // 0-1 representation of the matrix 

int best_ordering[MAX_SIZE];

double ngen; // variable containing the number of "nodes" generated
int matrix_size; // Actual size of the matrix based on input
int band; 	// bandwidth to achieve





// define max if necessary
#ifndef max
#define max(A,B) ((A)>(B) ? (A) : (B))
#endif

void print_ordering(int *);

int main(int argc, char *argv[]) {
	void minimize_band(void);

	int read_input(char *);
	int check_matex_symm(void);

	if (argc != 2) {
		printf("USAGE:   %s matrixfile\n\n", argv[0]);
		exit(1);
	}

	matrix_size = read_input(argv[1]);
	// check input 
	if (matrix_size == 0 || matrix_size > MAX_SIZE) {
		fprintf(stderr, "Illegal matrix file. Recall MAX_SIZE=%d\n", MAX_SIZE);
		exit(1);
	}
	if (check_matex_symm() == 0) {
		fprintf(stderr, "Matrix must be symmetric!\n");
		exit(1);
	}

	minimize_band();
	return 0;
}

/* ******************************************************************
 This function performs test in the following way: it builds
 a matrix in which the elements within the band.
 the matrix is then permuted
 randomly and bandwidth minimized using bida()
 ****************************************************************** */
void minimize_band(void) {
	void build_matco(void);
	void print_matex(void);
	void ida(result *);
	result ris;

	build_matco();
#if PRINT_MATRIX
	print_matex();
#endif

	ida(&ris);
	printf("Minimum bandwidth: %d\n", ris.band);
	printf("Time: %10.4lf(ms)\n", ris.elapsed_time);

}

/* ***************************************************************
 using Iterative Deepening Algorithm
 *************************************************************** */
void ida(result *ris) {
	int set_size(list);
	int reorder(int);
	int get_band(int *);
	int newb, bandwidth, i;

	//struct timeval t1, t2; // POSIX
	//time_t t1, t2; // Windows, get time now 
	clock_t t1, t2;
	double elapsedTime;

	// CLOCK START 
	ngen = 0;
//	gettimeofday(&t1, NULL); // Only on POSIX
    //t1 = _ftime64(0);
    t1 = clock();
	bandwidth = 0;
	for (i = 0; i < matrix_size; i++)
		bandwidth = max(bandwidth,(1+set_size(matco[i]))/2);

	for (; bandwidth < matrix_size; bandwidth++) {
		//gettimeofday(&t2, NULL);
		//t2 = _ftime64(0); 
		t2 = clock();
        //elapsedTime = (t2.tv_usec - t1.tv_usec) / 1000.0;
        //elapsedTime = difftime(t1,t2);
        elapsedTime = ((float)(t2-t1))/CLK_TCK;
#if USE_TIME_LIMIT
		if(elapsedTime >= TIME_LIMIT_AMOUNT) {
			printf("Time limit exceeded! \n");
			exit(1);
		}
#endif

		//printf("\nTry band:%d\n", bandwidth);
		if (reorder(bandwidth)) {
			//printf("FOUND!");
			break;
		}
	}
	// stop timer
	//gettimeofday(&t2, NULL);
	//t2 = _ftime64(0);
	t2 = clock();
	// compute and print the elapsed time in millisec
	//elapsedTime = (t2.tv_sec - t1.tv_sec) * 1000.0; // sec to ms
	//elapsedTime += (t2.tv_usec - t1.tv_usec) / 1000.0; // us to ms
    
    //elapsedTime = difftime(t1,t2); // for _ftime
    elapsedTime = ((float)(t2-t1))/CLK_TCK;
	
    
    printf("\nOptimal ordering:\n");
	print_ordering(best_ordering);
	newb = get_band(best_ordering);
	if (bandwidth != newb) {
		printf("Band verification failed! (%d vs %d) -ida-\n", bandwidth, newb);
		exit(1);
	}

	//  printf("CLK_TCK: %lf",CLK_TCK);
	//  ris->elapsed_time =  ((double) (end-start))/CLK_TCK;
	ris->elapsed_time = elapsedTime;
//	ris->node_gen = ngen;
	ris->band = bandwidth;

}

/* ****************************************************************
 This routine tries an ordering with bandwidth b.
 b = Upper(k/2), where k is maximum degree of the graph.
 **************************************************************** */
int reorder(int b) {
	int try(status *, Node, list);
	int try1(list, list, Node);
	list set_empty(void);
	int set_size(list);
	list set_add(Node, list);

	status start;
	int i, j;

	//initializes the global variables
	band = b;
	for (j = 0; j < matrix_size; j++) {
		start.queue[j] = j;
		start.expire[j] = band;
	}
	start.qsize = matrix_size;
	start.busy = 0;
	start.notfree = set_empty();

	//consider a node at a time
	int success = 1;
	list unselected;
	unselected = set_empty();
	for (j = 0; j < matrix_size; j++) {
		unselected = set_add(j, unselected);
	}

	for (i = 0; i < matrix_size - 1; i++) { // takes as the first node
		//printf("\nPermutation starts with: %d", i + 1);
		ngen++;
		if (set_size(matco[i]) > band)
			continue;

		//    if(try1(unselected, matco[i], i))
		//      return 1;
		if (try(&start, i, matco[i]))
			return 1;
		//    if(success==1) return 1;
		//    else continue;
	}
	return 0;
}

int try1(list unselected, list adj_list, Node n) {

	list set_difference(list, list);
	list set_union(list, list);
	list set_add(Node, list);
	list set_empty(void);
	int set_size(list);
	int set_member(Node, list);
	list set_intersect(list, list);
//	int check_compatible(list);
	if (n >= MAX_SIZE) {
		printf("OVER MAX DIMENSION -- try1");
		exit(1);
	}
	int j = 0;
	for (j = 0; j < matrix_size; j++) {
		int node = n + j;
		int k;
		list unselected1;
		list selected = set_add(n, set_empty());

		unselected1 = set_difference(unselected, selected);
		list adj_list1 = set_union(adj_list, matco[j]);
		adj_list1 = set_intersect(adj_list1, unselected1); // remove duplicated.
//		if (check_compatible(adj_list1)) {
//			if (try1(unselected1, adj_list1, n + 1)) {
//				return 1;
//			}
//		}
	}
	return 0;
}


/** Verifying compatibility
    if updated the number of busy nodes exceeds the bandwidth,
    the permutation will not make it. 
*/
int check_compatible(list adjlist, list toadd, status * new, int i) {
	list set_difference(list, list);
	
	Node nextnode = new->queue[i];

	adjlist = matco[nextnode];
	// Remove Duplicated, leave only diff between previous selection.
	toadd = set_difference(adjlist, new->notfree);

	// If current index reaches to the busy,
	// - all busy nodes are processed, then do not subtract.
	int newbusy = new->busy + set_size(toadd) - (i < new->busy ? 1 : 0);
	//printf("\nnewbusy %d", newbusy);
	if (newbusy <= band) return 1;
	else return 0;
}

/* ***********************************************************************

 By DFS, the routine will maintain the queue of nodes sorted by ascending
 expiration time order.

 When it picks a node, a node's adjacent nodes will be 'hold' (=busy, not free)
 while the exp time of them would be decremented.
 
 Less exp value nodes (busy) will be chosen first. 

  *********************************************************************** */
int try(status *old, Node n, list toadd) {
	void create_newqueue(status *, status *, Node, list);
	int init_status(status *old, status *new, Node n, list toadd);
	int try(status *, Node, list);
	int check_compatible(list adj_list, list toadd, status * new, int i);

	list set_difference(list, list);
	list set_union(list, list);
	list set_add(Node, list);
	int set_size(list);
	int set_member(Node, list);

	int maxn, i;

	status new; // the variable containing the next status
	Node nextnode;
	list adjlist;
	
	int newbusy;

//	printf("NODE %d", n + 1);
	// n.qsize initializes and checks to see if it is the last node
	new.qsize = old->qsize - 1;
	if (new.qsize == 0) {
		best_ordering[matrix_size - new.qsize - 1] = n;
		return 1;
	}

	// initializes the new status
	maxn = init_status(old, &new, n, toadd);
	

	// looking for a successor
  	for(i=0;i<maxn;i++) {
	    ngen++;                          // new node position generated
	    nextnode = new.queue[i];         // possible successor 
//	    printf("\nNEXT NODE %d %d/%d\n",nextnode+1,i,maxn);

	    adjlist = matco[nextnode];
	    // Remove Duplicated, leave only diff between previous selection. 
	    toadd = set_difference(adjlist, new.notfree);	
    
	    // If current index reaches to the busy,
	    // - all busy nodes are processed, then do not subtract.

	    int adjust;
	    if (i<new.busy) adjust = 1;
	    else adjust = 0;
	    newbusy = new.busy + set_size(toadd) - adjust;
	    
	    if(newbusy<=band) {	// After all busy nodes are processed, 
	      // passes to the next level
		if(try(&new, nextnode, toadd)) {
		        best_ordering[matrix_size - new.qsize-1]=n;
		        return 1;               
		}
	    }
	}     
	return 0;
}

/* *************************************************************
 initializes the new status on the description of the
 the old and the node and that was chosen. returns the number of
 nodes in new-> queue [] to consider
 ************************************************************* */
int init_status(status *old, status *new, Node n, list toadd) {
	void create_newqueue(status *, status *, Node, list);
	int get_max_index(status *);
	list set_union(list, list);
	list set_add(Node, list);
	int i, j, maxn;

	// initializes the list of nodes which are free
	new->notfree = set_union(old->notfree, toadd);
	new->notfree = set_add(n, new->notfree);

	// initializes the queue list []
	// this will generate new list which separates nonfree and free nodes
	// index will be new->busy (num of nonfree nodes)
	create_newqueue(old, new, n, toadd);

	// initializes the list exp_time []
	for (i = 0; i < new->busy; i++) { // for every busy nodes
		j = new->queue[i];
		// reduce each busy node's expire value by one. - for only nonfree(=busy).
		new->expire[j] = old->expire[j] - 1; // approaching expiration busy nodes
	}
	for (; i < new->qsize; i++)
		new->expire[new->queue[i]] = band; // free nodes have maturity = band
//	printf("node: %d", n + 1);
//	printf("\nqueue:\t");
//	for (i = 0; i < matrix_size; i++)
//		printf(" %d", old->queue[i] + 1);
//	printf("\nexp:\t");
//	for (i = 0; i < matrix_size; i++)
//		printf(" %d", old->expire[old->queue[i]]);
//	printf("\nqueue:\t");
//	for (i = 0; i < matrix_size; i++)
//		printf(" %d", new->queue[i] + 1);
//	printf("\nexp:\t");
//	for (i = 0; i < matrix_size; i++)
//		printf(" %d", new->expire[new->queue[i]]);
//	printf("\n");
	// calculates the number of nodes to be considered
	maxn = get_max_index(new);
	if (maxn < 0)
		maxn = new->qsize; // considers all nodes in queue []

	return maxn;
}

const char *byte_to_binary(int x) {
	static char b[9];
	b[0] = '\0';

	int z;
	for (z = 128; z > 0; z >>= 1) {
		strcat(b, ((x & z) == z) ? "1" : "0");
	}

	return b;
}

/* *********************************************************************
 This function takes as input a status (list of nodes busy and
 their expiration dates) and calculate how many nodes must be considered as busy.
 The function returns -1 if the nodes must also be considered not busy.
 Note: an expiration value of 0 indicates that the node must be chosen to
 this step. The nodes not busy have expiration == band because no adjacent node 
 is placed.
 ********************************************************************** */
int get_max_index(status *n) {
	int level;
	int i, counter;
	int next_time;

	// Printing nodes to consider
	//  printf("\nnodes to consider:\t");
	//  for(i=0;i<n->busy;i++)
	//    printf("[%d %d] ",n->queue[i]+1,n->expire[n->queue[i]]);
	//  printf("\n");
	

	counter = 1; /* counter, maximum number of nodes at this level */
	level = 0; /* initial level */

	for (i = 0; i < n->busy; i++) { // for every busy node
		next_time = n->expire[n->queue[i]]; //it's exp value
		if (next_time == level) { // exp==level?
			counter--; /* found new node at this level */
		} else { /* starts new level */
			if (counter == 0)
				return i; /* consider only the first nodes */
			else if (counter > 0) {
				counter += next_time - level - 1;/* num. nodes for this level */
				level = next_time; /* add level */
			} else {
				printf("Expirations inconsistent -get_max_index-\n");
				exit(1);
			}
		}
	}
	if (counter < 0) {
		printf("Expirations inconsistent -get_max_index-\n");
		exit(1);
	}
	if (counter == 0)
		return n->busy; /* considers only the nodes busy */
	else
		return -1; /* also considers non-busy */
}

/* *********************************************************************
 This function takes as input the old state and creates a new list
 queue. Variables are initialized new-> queue [], new-> qsize,
 new-> busy.
 New nodes are busy busy + the old ones toadd.
 ********************************************************************* */
void create_newqueue(status *old, status *new, Node n, list toadd) {
	int set_member(Node, list);
	int index_busy;
	int index_free;
	int i;
	Node qi;

	new->qsize = old->qsize - 1;
	index_busy = 0; // -- to occupy the first position in new->queue[] --
	index_free = new->qsize - 1; // -- last position to be filled in new->queue[] -

//	printf("\ncreate new queue\n");
	// init queue
	for (i = 0; i < matrix_size; i++) {
		new->queue[i] = 0;
	}
	// ------- copy nodes previously busy ----
	for (i = 0; i < old->busy; i++) {
//		printf(" %d", old->queue[i] + 1);
		if (old->queue[i] != n)
			new->queue[index_busy++] = old->queue[i];
	}
	//  printf("\n%s",byte_to_binary(toadd.bits[0]));
	//  printf("%s\n",byte_to_binary(toadd.bits[1]));
	//  printf(" / ");
	// ------- adds the nodes of busy old->queue[] -----------
	for (; i < old->qsize; i++) {
//		printf(" %d", old->queue[i] + 1);
		qi = old->queue[i];
		if (set_member(qi, toadd)) {
//			printf("busy[%d]\n", index_busy);
			new->queue[index_busy++] = qi; // puts qi in the busy (from front)
		} else if (qi != n) {
//			printf("free[%d]\n", index_free);
			new->queue[index_free--] = qi; // puts qi in the non busy (from the end)
		}
	}
	new->busy = index_busy;


}

/* ***********************************************************************
 performs the conversion from expanded format of the matrix to form
 compact
 *********************************************************************** */
void build_matco() {
	list set_add(Node, list);
	list set_empty(void);

	int i, j;

	for (i = 0; i < matrix_size; i++) {
		matco[i] = set_empty();
		for (j = 0; j < matrix_size; j++) {
			if (matex[i][j] != 0)
				matco[i] = set_add(j, matco[i]);
		}
	}
}

void print_matex() {
	int i, j;

	for (i = 0; i < matrix_size; i++) {
		for (j = 0; j < matrix_size; j++) {
			printf("%d ", matex[i][j]);
		}
		printf("\n");
	}
}

// return 1 if matex is symmetric, 0 otherwise
// also kill diagonal ones
int check_matex_symm() {
	int i, j;

	for (i = 0; i < matrix_size; i++) {
		matex[i][i] = 0;
		for (j = i + 1; j < matrix_size; j++) {
			//printf("%d,%d: %d \t %d,%d: %d",i,j,matex[i][j],j,i,matex[j][i]);
			if (matex[i][j] != matex[j][i])
				return 0;
		}
	}
	return 1;
}

/* ***********************************************************************
 calculates the bandwidth of the matrix matex [] [], when
 rows and columns are reordered with permutation p.
 *********************************************************************** */
int get_band(int p[]) {
	int i, j, pi, pj, d, band;

	band = 0;
	for (i = 0; i < matrix_size; i++) {
		pi = p[i];
		for (j = i + 1; j < matrix_size; j++) {
			pj = p[j];
			if (matex[pi][pj] != 0) {
				d = j - i;
				band = max(band,d);
			}
		}
	}
	return band;
}

/* *********************************************
 print permutation
 ********************************************* */
void print_ordering(int p[]) {
	int i;

	printf("{");
	for (i = 0; i < matrix_size; i++)
		printf(" %d", p[i] + 1); // starts from 0
	printf("}\n");
}

int read_input(char *infile_name) {
	FILE *file;
	int k, nodes, edges, i, j;
	if (infile_name == NULL)
		return 0;

	file = fopen(infile_name, "r");
	if (file == NULL) {
		perror("Unable to open matrix file!");
		return 0;
	}

	fscanf(file, "%d", &nodes);
	printf("nodes: %d\n", nodes);
	fscanf(file, "%d", &edges);
	printf("edges: %d\n", edges);

	if (nodes > MAX_SIZE) {
		fclose(file);
		return 0;
	}

	while (!feof(file)) {
		fscanf(file, "%d", &i);
		fscanf(file, "%d", &j);
		matex[i - 1][j - 1] = 1; //start from 0
		matex[j - 1][i - 1] = 1;
	}

	fclose(file);
	return nodes;
}


/* ===============================================================
 For bit manipulation
 These assume that list is an array of a structure contentnente lSize
 unsigned long (an element I 'identified by an integer
 between 0 and lSize-1 (== n/32), and an integer between 0 and 31 (n == 32%))
 =============================================================== */

list set_intersect(list a, list b) {
	int i;
	for (i = lsize - 1; i >= 0; i--)
		a.bits[i] = a.bits[i] & b.bits[i];

	return a;
}

list set_union(list a, list b) {
	int i;

	for (i = lsize - 1; i >= 0; i--)
		a.bits[i] = a.bits[i] | b.bits[i];

	return a;

}

int set_member(Node n, list a) {
	int q, r;

	q = n / 32;
	r = n % 32;

	if ((1L << r) & a.bits[q])
		return 1;
	else
		return 0;
}

list set_add(Node n, list a) {
	int q, r;

	q = n / 32;
	r = n % 32;
	a.bits[q] = a.bits[q] | (1L << r);
	return a;
}

list set_empty() {
	list a;
	int i;

	for (i = lsize - 1; i >= 0; i--)
		a.bits[i] = 0L;
	return a;
}

list set_difference(list a, list b) {
	int j;
	unsigned long i;

	for (j = lsize - 1; j >= 0; j--) {
		i = a.bits[j] & b.bits[j];
		a.bits[j] = a.bits[j] ^ i;
	}
	return a;
}

int set_size(list a) {
	int set_size32(unsigned long);
	unsigned long u;
	int i, sum;

	sum = 0;
	for (i = lsize - 1; i >= 0; i--) {
		if ((u = a.bits[i]) != 0)
			sum += set_size32(u);
	}
	return sum;
}

int set_size32(unsigned long a) {
	static int first_call = 1;
	static char ones[256];
	int i, j, a_size;

	if (first_call) {
		first_call = 0;
		for (i = 0; i < 256; i++) {
			ones[i] = 0;
			for (j = i; j != 0; j = j >> 1)
				if (j & 1)
					ones[i]++;
		}
	}

	a_size = 0;
	while (a != 0) {
		a_size += ones[(int) a & 0xFF];
		a = a >> 8;
	}
	return a_size;
}
